# -*- coding: utf-8 -*-

"""
    2025-02-01 by Jie Zheng
"""


import numpy as np


def match_radec(ra1, dec1, ra2, dec2, dislimit=0.001):
    """
    match two sets of ra,dec, transfer them to xsi,eta, and match them
    they should be in a small area
    :param ra1: ra1, in degree
    :param dec1: dec1, in degree
    :param ra2: ra2, in degree
    :param dec2: dec2, in degree
    :param dislimit: limit of distance, in degree
    :param crval: crval, in degree, if None, use the mean of ra1&2,dec1&2
    :return: ix1 and ix2, the index of matched pairs
    """
    # number of points
    n1 = len(ra1)
    n2 = len(ra2)
    # transfer degree to radian
    ra1r  = np.deg2rad(ra1)
    dec1r = np.deg2rad(dec1)
    ra2r  = np.deg2rad(ra2)
    dec2r = np.deg2rad(dec2)
    dislimit_line = 2 * np.sin(np.deg2rad(dislimit / 2))
    # ra,dec to xyz (assume r = 1)
    x1 = np.sin(dec1r) * np.cos(ra1r)
    y1 = np.sin(dec1r) * np.sin(ra1r)
    z1 = np.cos(dec1r)
    x2 = np.sin(dec2r) * np.cos(ra2r)
    y2 = np.sin(dec2r) * np.sin(ra2r)
    z2 = np.cos(dec2r)
    # distance between two point sets
    # distance matrix of each x1-x2 and y1-y2
    dx = np.array([x1] * n2).T - np.array([x2] * n1)
    dy = np.array([y1] * n2).T - np.array([y2] * n1)
    dz = np.array([z1] * n2).T - np.array([z2] * n1)
    # 2d distance (squared)
    dis = np.sqrt(dx * dx + dy * dy + dz * dz)
    # find pairs close enough
    ix1, ix2 = np.where(dis < dislimit_line)

    return ix1, ix2

